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Functions for Thermal Stress Calculation Near a 
Transient Heat Source on a Flat Surface 
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When an initially unstressed elastic solid at uniform temperature is subject to a transient, locally 
two-dimensional heat flux on a flat surface, the two-dimensional total stress field near the wall is locally 
determined for short times and may be constructed from the functions described in this report; For- 
tran programs are available for their computation. In particular, maximum total stress and total stress 
fields for initially large heat fluxes are readily obtained for estimations of yield probability. 
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1. Introduction 



When heat is suddenly applied to or withdrawn from 
the surface of an elastic solid, large superficial stresses 
are set up which can cause surface damage when local 
elastic limits are exceeded. This occurs, for example, 
when hot or cold material is suddenly brought into 
contact with the solid surface. If the surface is reason- 
ably flat and the initial heat flux distribution rapidly 
varying in one direction only (for example, near the 
perimeter of contact of the hot material mentioned 
above), the initial transient stresses near the wall are 
obtainable from plane strain theory. 

In this report, plane strain theory is used to analyze 
the total stress field near the surface of a semi-infinite 
elastic solid during a short time interval during which 
a large two-dimensional heat flow occurs; formulae 
are presented for numerical evaluation of the stress 
for a class of heat flux pulses which may be approxi- 
mated by a simple physically relevant form. Such a 
stress field is useful for estimating the probability 
and location of fracture from either maximum stress- 
yield strength data or a more sophisticated probability 
theory. 

We shall consider then a semi-infinite elastic solid 
occupying the region y^O which at time £ = has 
vanishing stress components, (7y = 0, and uniform 
temperature T=0. Cartesian coordinates x and z are 
taken to lie in the surface and a given outward heat flux 
q(x, t) independent of z is stipulated, so that the tem- 
perature T{x, y, t) , the stresses (Tij(x, y, t) , the strains 
€ij{x, y, t), and the displacements m{x, y, t), (i = l, 2), 
do not depend on z. We assume also (plane strain) 
u,i(x* y, t) = 0. 



We shall neglect the heat generated in the solid 
by the strains, and assume that the temperature is 
governed by the heat equation 



r,,= Z)(7\,, r +7\ ,,,,), 



(1) 



where D is the coefficient of thermal diffusion. We 
shall also neglect the inertia] effects in the solid, as- 
suming that the stress is quasi-stationary near the 
surface. Then 

o-ijj = (), (i,7=l,2,3), (2) 

where the stress-stain relation [1] 3 is given by 

dij = kdij€ kk H- 2fxeij — rriTSij (3) 

the strains are given in terms of the displacements by: 

£u = h( u iJ+ u J. 0, (i,/=l, 2, 3). (4) 

Here k and fx are the Lame constants, a is the coeffi- 
cient of linear thermal expansion, and m = a{3k + 2/x). 
The Lame constants k and /jl are given in terms of the 
Young's modulus E and Poisson's ratio v as follows [ 1 ] : 



vE 



E 



(i + ^)d-2^r 



m : 



2(l + i/) 



The validity of the assumptions made in the above 
equations will depend on the specific heat flux in- 
volved, q(x 9 t), and the region being studied, and should 
be checked a posteriori. 

We shall assume further that a single scalar, the 
total stress (the first stress invariant), 



O" = (Tu = <J\ i + 0~22 ~r~ 0"33, 



(5) 
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suffices for yield analysis. Yield criteria may be The equations (2) become: 
based on extreme values of a or on the more sophis- 
ticated probability theories using the stress distri- {k-\-2fx)ui, xx + ku 2 , xy — mT, x + f*<(ui, yv + u 2 , xy ) = 
bution [2, 3]. , /ir>x 

Boundary conditions on the stress system are that 

surface forces vanish on y^O: (X + 2/x)w 2 , yy + kui } xy — mT, y -\- fi(ui, xy + ii2, xx) = 0. 

' ' For a function f(x, y, t), satisfying appropriate con- 

T+ . ii t r A -\ .i . i . 11. . i-rr . ditions, we introduce the Laplace transform Jf in t 

It is well known 4 that when two bodies at different , , r c ^Y 

L J . . . . , nu *■* an d the rouner transform & in x: 
temperatures are put in contact at £ = U, heat is trans- 
ferred at a rate 0(£~ 1/2 ) for sufficiently short times. F( . l. \ = <& o?ft t \ 
We shall assume the heat flux time history in a time ^' ' ' J\ -> y-> ) 



(13) 



interval of interest, ^ t ^ £f, is approximated as 1 f* f«= 

follows: = ^T dxe dte-«f(x, y, t) 

3 

q(x, 0=2 ^nM* n_1/2 , (0 ^ t ^ t F ). (7) with the inverse: 

n = 

In particular, we shall find the total stress fields 
(Tndix, y, t) associated with the line source heat fluxes: L_ I * dke Xkx 

q n8 (x, t) = d(x)t»-*l\ (* = 0, 1, 2, 3). (8) 

fc+/oc 

Then the stress <r(x, y, t) for the general heat flux of 2ni J c -too 

(7) is given by 

J^ f x The transforms of (u u u 2 , T, a) are denoted, respec- 

cr(x, y, 0= 2 J dsA„ILs)<rndix — s* J, *)■ (9) tively, by (£/, F, 0, 2). 

" =0 x We obtain the system of ordinary differential 

equations: 



Our principal task is to construct the functions 
o-ndx, y, t). 



2. Analytic Evaluation of the Stress (Tndx, y, £) 

For the case of plane strain, the strains and stresses 
are given in terms of the displacements as follows: 

€n =1*1, x 

ei2 = € 2 l = i(Ui, y + U2, x) 

£22 = y>2. y 

e i3 = e 3i = (j=l, 2, 3) 



(s + k 2 D)Q- DQ, yy = 

-(k + 2fjL)k 2 U + nU, yy + (\ + fjL)ikV, y -ikmQ = 0, 

(15) 
{- tik 2 V+(k + 2iJL)V, uy + (k + rfikU, y-mG, y = 0. 



Boundary conditions are boundedness of the trans- 
forms as y— » oo, and the following conditions at y=0 
follow from (6): 



U, y +ikV=0, ikkU+(k + 2/jL)V, y -mG = 0, (16.1) 
( 10 ) and by (8) at y=0: 



eu = u U x + u 2 , y °' y= ^~^ (16 - 2) 

o~i i = ken ■+" 2/xwi , j- — ml 

o~i2 — cr 2 i =ti(Ui, y -\- Uz, x ) where K is the coefficient of heat conduction. The 

system of (15) has solutions of the form 



o "22 — ken + 2/ULU2, y — mT 
0-33 = keu — mT 
0-23 = cr 32 = or* = 0-31 = 
cr = (Tu = (3k + 2jLt)€// — 3/717 1 . 




(ID /e\ /eA /o 

f/, |e .«+[ f/ 2 +y^3 I 
Vj \V 2 +yV 3 , 
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V\ 



k 2 + 



I) 



1/2- 



e— sgn A, 



(18) 



Inverse transforms of the terms in (21) may be ob- 
tained by manipulation of tabulated transforms. It 
is convenient to introduce the variables: 

.(gj y» r) 



and v\ has nonnegative real part. Substitution of 
(17) into (15) and use of the conditions (16.1) gives 
the following: 



(p = tan 



r=vx 2 -\-y 2 , (f, 17, p) : 



T? 



z = p 2 e~ 2lV , 



/4D* 



arg z ^77. 



(24) 



F f/=biie-^-h(gi2+ygi 3 )e- u ^]ei 
. F = [fti«-*rf+ (gr 22 + yg 23 )e- ,w »]e 1 , 



(19) 



where: 



<7n : 



421 : 



912-^7 



ikmD 
"(\ + 2/i,y 

v\inD 
mD 



Then for the first of the terms in (21), we note [5a]: 

L = j«-i/ % -w^/l +lli I- | z |\ (25) 

27T V^ / 



where we have made use of the relation between the 
Whittaker function W- n .dz) and the confluent hyper- 
geometric [6a] function i//(* +n, 1 ; 2). 



i^- 1 ^- 



eft+3/2 



1 



f n-l/2 



cos 2<£> 



V2^ r / + g H 



(26) 



<?13 r 



<?22 = 



2iks(k + n)(\ + 2fi) 

iml) 
'2/jls(\ + 2ijl) 

mD 



{2e(\ + 2fjL)kp l -2fik 2 }, 



(20) 



2&jMA +//,)( A + 2/x,) 

^23= J€^i 3 . 



{2/J^,-2^A- 2 e}. 

{- 2/x 2 A^i + 2//,(A + 2^)A 2 e} , 



For the third term of (21) we note that 



j, i5 /< + .}/2 Y^ ^ 



(27) 



^here: 



1 



/„(*, y, - ^~ x —^~ 2 J(x, y,s), 

</(*, y, s)= dk 
Jo 



s «-r3/2 

cos A;.% • e -A ^ 



(28) 



The stress transform £ is readily found to contain no 
term linear in y: 



k 2 + 



D 



2 = 



I /• I /-2' 



e, 



Then in terms of Struve's function H</z) and Bessel's 
function [5c] of the second kind KoU), 



(21) 



J(x, y, s)- 



where the a,- do not depend on (A, 5, y): 

4/xw 
X-h2/x 

2/Lt///.P(3\-r-2jLt) 
(X + /LL)(X + 2ft)' 



^"M^VI 



x /y-y (^^ 



(29) 



«:> = — CK3 = 



(22) 



We notice that the first term on the right-hand side of 
(21) gives a term directly proportional to the tempera- 
ture (17). By (16.2) and (17), 

1 Wn±h) 

e i= (23) 



Then [5d, 6b], 

/-i(«, y, *) = — 7= Real \e^K ft 



= ^Real{*(|,l;; 



(30) 



a:V27t ^,5" +i/2 



vhere Ko is the modified Bessel function of the second 
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kind. It follows that [5e]: 



/„(*, y, t) = ( J 



An(f,T?) = 



Hti- 



d* 



e-*X»+* f l;p 2 ). (35.2) 



/-i(x, y, *) 

^ \n+l/2 



; 2 ( ->" + ' 4D 1 Rea1 ' 



-ip(2n+l) 



/;r'^(i-)i 



(31) 



The constants A, in (34) are the following: 
/x(3X + 2 A t) 

A 2 : 



"2(A + /x)(A + 2p,)' 
2/x 



(36) 



(X + 2p,) 

where (fdt • )" +1 is the (ra-hl)-th iterated integral. 

Each integral is independent of the argument of the We note the useful relation: 
lower limit for z if |arg z\< 37r/2. It is easy to estab- 
lish by recursion the following formula [6e]: 



(L$r*&»)-&*(—i* 



r »+: 



n ( — \k r* 



The finite sum removes from the first term of the 

right-hand side those terms of the asymptotic series 

which would give a nonvanishing result as z— > °° [6d]. 

We introduce the following real- valued functions: 

Gm(£, 17) = Real wf^ - ' 1 ' l '* z )\* 
Gtn({, rj) = Real W - l)zi// (J- n 9 3; z) 

(33) 



+ (2n-hl)i//(--n, 2;z 



Finally then we are able to write the total stress as 
follows: 



crUx, y, t)=-t n -^{KJ n {t, 7l) + AJU(, 71)} (34) 



w}iere the functions (f n , h n ) are even in £ and do not 
depend on the material parameters: 



r/> v 12 COS 2<f 



(-)" 



(2n+l)r U + 



n (_)K 
■\Gtn(£, tj)~2 "jFT 



(2n+ 1 - 2K)(2n - 2K)p 2n ~ 1 - 2K cos p (2n- 1 - 2A) k 



r n+ 



r [n-K+\ 



T„ s (x, y, i)-- 



J_ 
2K 



t-^Wf , rj), 



(37) 



giving the temperature field in the solid due to the line 
heat source of (8). 



3. Numerical Evaluation of (/„, hn) 

Expressions are given in appendix A for the interim 
functions Gi n (£, 17), G2n(£, 17), both convergent and 
asymptotic series. Using double precision computa- 
tion (25 digit mantissa) for the convergent series, 
single precision computation (10 digit mantissa, 
exponent less than 308) for the asymptotic series, 
agreement between the two was to six significant 

figures at p 2 = 20 for several values of <p on ^ (p ^ — , 

and thus p = V20 was taken as the crossover point. 

Expressions for /i n (£, 7}) are also given in appendix 
A in terms of the functions Gm(0, p) using a recursion 
relation. 

As p - »0, the G2n have integrable singularities: 

/, -> J {in p 2 + (^ (I- nj - i#l) - #2) + 2) + 0(p)}, 
^^-i{lnp 2 +(i//(|+/i) -2i/<l)) +0(p)}, 



(38) 



where 0(p) vanishes as p — » 0, and i//M is the Digamma 
function [7] . We may then integrate a n s over y = as 
required in the convolution (9). 
For large p, the simple forms, 



/. — ± 



1 2 cos 2<p 



r(»+ 



tt2« + 1 



77T(n + 2) p 



i2253£+0(^, 



A. ~ 0, 



(39) 



J give better than one percent accuracy for p ^ 10. 
(35.1) These formulae show that 
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Vnb 




t n y 

2 + y 2' 



(40) 



so that for short times (r/v^Dt —> °°), the stress field 
of a line source dies off 0(l/.x 2 ) for \x\ >y (although 
only 0(l/y) for y> \x\). This integrable decay rate 
in x makes the convolution (9) useful since the large 
stresses expected near the surface at early times are 
locally induced, and the plane stress theory is applica- 
ble in more general situations. 

A subroutine for the calculation of {Gm(^, 77), 
G2n(^, T?)}, and a program employing this subroutine 
for the calculation of {/ w (£, rj), h n (^, 77)} for a given 
value of 77 on the set of values, 



£ = 0.0(0.1)10.0, 



(41) 



_is_dis cussed in appendix B. The calculated values 
are stored on tape in consecutive values of £ in the 
following order for each f : 



{f (i,v), ■ ■ .,/•(£, i|),M6i) 



., M£, v) 



(42) 



The values are sufficiently slowly varying that Simp- 
son's rule integration in (9) is normally adequate. 
Calculation of these functions on the set (41) for a given 
value of 1) takes on the average about two minutes, 
and was done for the set of integer values 



77 = 0.0(1.0)10.0. 



(43) 



These data are now on magnetic tape, and are available 
on punched cards. (The singular values (£ = 0, 17 = 0) 
have been replaced by the values at (£ = 0.1, 77 = 0) as 
being adequate for integration.) 

4. Step-Input Stress Response 

A special case of heat flux of particular interest for 
stress evaluation near the perimeter of an (almost) 
uniform source is obtained from the general source of 
(7) by setting 



A n (x)=A„h(x), 
where the A n are constants and 

f (x<0) 

h(x) = 

{ 1 (*2*0). 

Then the resulting stress cr/, is found by (9) 



(44) 



(45) 



to be 



<r h (x 9 y 9 t)=- 



K 



1 = 

+ \,H n (f.n)}, 



(46) 



yhere: 




(M\s\ 


u) 


\hn(\s\ 


,v) 


(x,y) 




VWt 





(47) 



(£ v) ■ 



When the entire path of integration (£1, £ 2 ) lies out- 
side the circle of radius p=10, the asymptotic forms 
(39) are used to evaluate: 



<M6,6;i>) : 



■\: 



dsfn(\s\, V ), 



each term of (39) having a simple closed-form 
Then for ^ 77 ^ 10, we have 



(48) 



ntegral. 



10.: 




<M£>~°°; y) 







(49) 



10.: 




(50) 




F„(10,ij)> 

Hn(l0, T?) i 



+ 



frtf, 10; v )\ 

} 



(51) 



for 77 > 10, (49) is valid for all f. 

Appendix C describes a program for determining 
{/%„ H „} for given 77, (one of the set (43) of course) on 
the set 



f= 10.(0.2)20., 



(52) 



using Simpson's rule for the evaluation of the integrals 
in (50). These results are on magnetic tape and are 
available on punched cards. Calculations took about 
one-half minute for each value of 77. 

When Aq=\, A l =A 2 = As = 0, the total stress field 
is self-similar in time, dependent £ and 77 only, and 
rises monotonically (dominated by G n (i;, 77)) in £ to a 
maximum value at £ = °° independent of 77. This value 
is attained at (77 = 0, £ = 10). to the accuracy of the 
calculation, and may be determined from the values 



fF (20, 0)\ 

Wo, o)J 



-0.0637^ 



1.7724 



(53) 
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For large values of £, the stress field (Th given by (46) 
is uniform in x for any coefficient set {A n }, and agree 
with the results for heat flux q(t) uniform in x described 
in appendix D. 

5. Conclusions 



7. Appendix A. The Functions {G ln {£;, 77), 

G2n(£, I?), h n {^ 7])} 

The following expansions for G ln , G 2 « are readily 
obtained from the definitions of (33) and the properties 
[6c] of the confluent hypergeometric functions i//(a, b; z): 



Formulae are presented for the estimation of total 
stress fields near the surface of a solid during a short 
time interval during which a large heat transfer occurs. 
Such a stress field is useful for estimating the proba- "in(g? V) z 
bility and location of fracture, either from maximum 
stress-yield strength data or from a more sophisticated 
probability theory. The stress fields obtained may 
be expected to be valid in a region adjacent to a point 
on the solid surface of dimension / ~ V^Dt and for the 
time interval ^ t ^ tp, the period of validity of the 
heat flux approximation (7), so long as the following 
orders hold: 



^Dt<{L, R m ,R q }, 



(54) 



■ n~\- r 



-1 \^ 1 ' \2 _ 

H^ (r!,! r (H) 



• p 2r [(ln p 1 ) cos 2np — 2ip sin 2r<p] 






iM|-7i + r)-20(l + r) 



where L is a length scale for the solid, and R m and 
R (} are, respectively, the mean radius of curvature of 
the solid surface and the radius of curvature of lines 
of constant outward heat flux at the point in question. 
It would be desirable to have a "maximum princi- 
ple" available, stating, for example, that if q(x, t) = 
for t > If, then the maximum value of cr(x, y, t) over 
(it, y) for any t > tF does not exceed the maximum value 
of <r(x, y, If) over (x, y). The part of the stress associ- 
ated with h n {£;, rj) snares with the temperature field 
given by (37) such a maximum principle, but no more 
general result appears available. 



• p' lr cos 2r<p\, (Al) 



\ 1 n 2n - 1 - 2m 

Gm(£,v)^ > H mB r— 

m + --n 



cos (p{2n— 1 — 2m) 



r {2~ 



TT 



IpI^tt/- 
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c 2 „(£, v)= { -^^ r (" + i) (2n+ 1} ln p 1 



^|-/7)-^(l)-^(2) + 2)+ 2 C °\ 2ip 
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irl r(|-njr(Z+2) 



[(In p-) cos 2l(p — 2(p sin 2lip] 



+ (1 + 



2n) > - 

4ril\ 



m-n + / 



jr? lW+2)T(~n 



2 + (2/+l) 



$Q-n+l\- 0U + /)-«K2 + /)) 



p-> cos 2Up , (A3) 
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Gue, v) 



00 

I 



ml 



(2n + 1 — 2m) (2n — 2 m) • p- 



' cos <p{2n — 1 — 2m), 



<P\ 



(A4) 



Let 



«(£,!))- 



r ln + 



e"C',„(0,p), 



(n=l,2). 



(A5) 



Then A»(£, 77) are found from the following expressions: 
ho=h$, 



hi=2p%* 



•4/if, 



h, = -[4p\2+p')-l]h* 



-|(2+p 2 )At, 



(A6) 



A 3 = 



10 



-<4 + p^(4p^(2 + p^)-])-4p- 



A* 
"0 



-(4 + p 2 )(2-fp^)-l 



/if. 



8. Appendix B. Description and Use of Fortran 
Programs for Calculation of {f u , h n ] 



vides the input values £ and 77 and writes the com- 
puted functions on a tape. Subroutine FSGSHS per- 
forms the actual computation of the functions (f n , h n ) 
for each value of (£, 77). Subroutine GSUMS com- 
putes the functions {Gm(£ , 77), G 2n (£; , 1?)} for each (£, 77) 
using the coefficients of G Ut and G 2w as computed in 
the program GCOEFFS. Subroutines GCOEF1 and 
GCOEF2 set up tables of these coefficients for use 
by GSUMS. 

The first input variable, ISKIP, is used to position 
the output tape, giving the number of files to be skipped 
before any writing is done on that tape. Since the 
functions {f n , h n } are computed for £ = -10.(0.1)10. 
for each 77, it is necessary to read in only the variable 
r](DETA). The program will read in an 77, do the 
computations and return to read the next 77, continuing 
until no cards remain. Note that 77 is read in as a 
double precision variable. 

Program listings for the CDC 3600 may be obtained 
from the authors. 

9. Appendix C. Description and Use of Fortran 
Programs for Calculation of {F n , H n } 

The program INTEGFGH computes the functions 
{F,„ H n } using the formulas (49), (50), (51). The 
functions {f n , h n } are read from the tape generated 
by program FORMFG. Results are written on an out- 
tape in the same order as are the {/„, h n }. 

Input to this program is in two parts: the tape gen- 
erated by FORMFG, and cards. The first card input 
variable is KSKIP, which positions the input tape by 
skipping KSKIP files. This option makes it possible 
to select the value of 77 for which the computations 
are to be performed. The upper limit of £(XIF) and 
an approximation to — o°(AA) are also card input pa- 
rameters. After each complete integration the input 
tape is rewound and the next value of KSKIP is read 
in, followed by the upper limit of £ and the approxi- 
mation to — 00. The program continues until no cards 
remain to be read. 

Program listings for the CDC 3600 may be obtained 
from the authors. 



Appendix D. Heat Flux Uniform in 



From the expansions for G\ n and G^n, given in ap- 
pendix A, it can be seen that there are multiplicative 
factors in the sums which are dependent only on n 
and the index of summation. These factors were 
therefore computed using the program GCOEFFS 
and were punched on cards in a format suitable for 
use directly as part of an object deck in another pro- 
gram. It was decided that the maximum number of 
terms necessary for required accuracy is 20 for the 
asymptotic expansions and 100 for the convergent 
series. 

Computation of the functions {/»(£, 77), h n (^, 77)} is 
done in program FORMFG using subroutines GCOEF1, 
GCOEF2, FSGSHS, and GSUMS written in CDC 
3600 Fortran. The main program FORMFG pro- 



When the heat flux of the wall in independent of x, 
the stress field will be also. Simple results may be 
obtained for this case, and these results have been 
used to corroborate the numerical results described 
in the text. 

Putting k = in (21) gives the ^-independent case 
for the Laplace transforms 



IT 



(Dl) 



If the outward heat flux q(t) has Laplace transform 
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Q(s), we have 



Q(s) = -K^G t (s), 



A — a,. 



(D3) 



so 



2 (* 5 > ; 



a.\e 



-K 



Q(s). 



A useful relation can be derived for the wall stress 
cr(0, t) in terms of the wall temperature T(0, t) only, 
independent of the heat flux history. Setting y=0 in 
(Dl) leads to the relation 



cr(0, i)=aiT(0, t). 



(D4) 



Using the convolution theorem for Laplace trans- The relations (D2) and (D4) are preserved by numeri- 

forms, we may obtain the following formula: cal tests of (46). 



a(y, t)=-2 x - 

17T 



-A dpq( 

7T Jo 



t(l-p 2 ))e-p 2 , (D2) 



(Paper 70B4-188) 
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